Seismic attenuation can be quantitatively parameterized by the quality factor Q, a parameter assumed to be constant over the frequency range used in geophysical applications. Accurate estimation of Q distribution and Q values is critical for geophysical exploration and production, such as rock property characterization, reservoir development, and attenuation compensation in earth imaging.
Q tomography is an advanced approach for automatic estimation of subsurface Q anomaly geometries and the associated Q values for 2D or 3D scenarios. This approach analyzes the attributes of the recorded seismic data to reconstruct the Q profile. Generally, Q tomography algorithms are classified into two main categories. One category is ray-based tomography (Quan and Harris, 1997; Rossi et al., 2007). The other category is wave-equation-based tomography (Liao and McMechan, 1996; Pratt et al., 2003). Wave-equation-based tomography is physically more accurate but computationally expensive and not practical for 3D cases. The present invention belongs to the category of ray-based Q tomography.
Ray-based Q tomography is essentially a linear optimization problem. Three main components of ray-based Q tomography algorithms are: 1) construction of a kernel matrix by using the raypath information obtained through ray tracing procedure; 2) construction of a measurement vector by using one or several attributes of the recorded seismic traces because these traces carry a great amount of information of the subsurface Q distribution; 3) solving the linear optimization problem formulated by relating the kernel matrix, the Q distribution profile, and the measurement vector. Most existing ray-based Q tomography algorithms use transmission seismic data to reconstruct subsurface Q models (Quan and Harris, 1997; Hu et al., 2011). This type of Q tomography algorithms are referred to as refraction seismic data Q tomography or transmission seismic data Q tomography. The procedure of transmission seismic data Q tomography is relatively simple because pre-migration pre-stack seismic data (i.e., time domain data) are used (Hu et al., 2011). Consequently, the construction of the kernel matrix for transmission seismic data Q tomography is straightforward. First, with the given seismic survey geometry and the given velocity model, shot-to-receiver ray tracing is implemented. Then, the number of seismic rays penetrating each subsurface grid is populated and the lengths of these rays in each grid are measured. Collecting the raypath information, one can build a kernel matrix. Combining the kernel matrix and the measurement vector, eventually one can formulate an optimization problem to reconstruct the subsurface Q model.
Unfortunately, transmission seismic data are not always available. In many exploration geophysical applications, the number of transmission seismic traces is very limited, which implies that the formulated transmission seismic Q tomographic inverse problem can be very underdetermined. Furthermore, with limited range of offsets, transmission seismic rays do not travel down into deep regions. In other words, these rays only carry the shallow region geophysical property information. Therefore, the transmission seismic data cannot be used to reconstruct deep Q models. FIGS. 1A-1B give an example of transmission seismic ray coverage and its relationship to the Q models at different depths. FIG. 1A shows the velocity model and the corresponding transmission seismic raypaths. In FIG. 1B, there are two Q anomalies. One of them is located in the shallow region, while the other is buried in the deep region. The shallow Q anomaly is fully covered by the transmission seismic raypaths. Therefore, this Q anomaly can be reconstructed by using transmission seismic data Q tomography. However, the other Q anomaly is too deep to be recovered by using transmission seismic data only. On the other hand, reflection seismic rays do penetrate this deep Q anomaly, as shown in FIG. 1B. Therefore, to estimate the deep region Q distribution profile, reflection seismic data have to be utilized. This invention is a method for using reflection seismic data in Q tomography to reconstruct a subsurface Q model reliably, especially for deep regions.
In transmission seismic data Q tomography algorithms, source-to-receiver ray tracing is performed and then the raypath information in each model grid is obtained to build the kernel matrix for Q tomography. The reason that this procedure can be performed in transmission seismic data Q tomography is that the raypaths are relatively simple for transmission seismic data, as shown in FIG. 2. It is straightforward to measure ray lengths lij in each grid to build the kernel matrix; for example l13, l14, l24, l25, l26, l36, l37 . . . in FIG. 2. In reflection seismic data Q tomography, usually the kernel matrix cannot be constructed by directly measuring the seismic raypath for each pre-migration seismic trace due to the following reasons: 1) there are too many reflection seismic traces; 2) reflection seismic raypaths are too complicated and it is extremely difficult to trace all the shot-to-receiver seismic rays; 3) there are many multiple arrivals which makes it impossible to separate the contributions from different seismic raypaths; 4) it is very difficult to pick the right events in the pre-migration reflection seismic traces for a specific raypath. For example, in FIG. 3, where the velocity is constant and only two reflectors are in presence, for the shot-receiver pair shown in FIG. 3, there are at least four different raypaths labeled as Ray 1 to Ray 4. For more realistic scenarios, the raypaths are even more complicated. Under some circumstances, if the geological structures are simple, it might be possible that pre-migration reflection seismic data can be analyzed and used for Q tomography (Rossi et al., 2007). However, in most cases, reflection seismic data Q tomography algorithms only work on post-migration seismic data (i.e., image domain data or depth domain data) instead of pre-migration seismic data (i.e., time domain data).
In most existing reflection seismic data Q tomography algorithms (Hung et al., 2008), the relationship between the post-migration seismic data and the subsurface Q distribution profile is established through spectral ratio method, which is a widely used approach utilizing the amplitude decay information to derive Q values. However, the amplitude decay based method for Q estimation is significantly affected by reflection and transmission loss and an illumination imbalance issue. Another method to link the seismic data and the Q profile is the so-called centroid frequency shift method (Quan and Harris, 1997), which is believed to be more robust since it is independent of geometrical spreading, reflection and transmission loss, and illumination imbalance. Unfortunately, the conventional centroid frequency shift method is only applicable to pre-migration seismic data (i.e., time domain seismic data). He et al. (2012) propose an approach to obtain correct spectral information from post-migration depth domain seismic data and then input the spectral information to the centroid frequency shift method. However, their method is applicable only to common angle gathers of seismic data, where the takeoff angles at the selected CDPs are known.